/*
 * Propeller.cpp
 *
 *  Created on: Apr 12, 2010
 *      Author: andrew
 */

#include "Propeller.h"
#include <cmath>
#include <algorithm>
using namespace std;


Propeller::Propeller(int Nb, double radius, double hub_radius, double root_chord, double tip_chord, double root_inc,
			double tip_inc, double cd_0, double cd_1, double cd_2, double cl_0, double cl_max){
	actualRadius = radius;
	Nblades = Nb;
	hubR = hub_radius;
	rootC = root_chord;
	tipC = tip_chord;
	rootI = root_inc;
	tipI = tip_inc;
	cd0 = cd_0;
	cd1 = cd_1;
	cd2 = cd_2;
	cl0 = cl_0;
	clmax = cl_max;
	J = 0;
	eta = 0;
	CT = 0;
	CQ = 0;
	CP = 0;
}

double Propeller::getRadius(){
	return actualRadius;
}

void Propeller::setAdvanceRatio(double Jnew){
	J = Jnew;
	// Dimensionless version
	double RPM = 1;
	double RPS = RPM/60;
	double radius = 1;
	double omega = RPM*2*M_PI/60;
	double U = J*RPS*2*radius;
	double dr = (radius-hubR)/Nstations;

	double r[Nstations], chord[Nstations], incidence[Nstations], du[Nstations];

	for (int i = 0; i < Nstations; i++){
	    r[i] = hubR + (radius-hubR)*(i+1)/Nstations - dr/2;
	    chord[i] = tipC*r[i]/radius + rootC*(radius-r[i])/radius;
	    incidence[i] = tipI*r[i]/radius + rootI*(radius-r[i])/radius;
	    du[i] = 0;
	}

	for (int it = 0; it < maxit; it++){
	    double T = 0;
	    double Q = 0;
	    for (int j = 0; j < Nstations; j++){
	        double V = sqrt(pow(U+du[j],2) + pow(omega*r[j],2));
	        double phi = atan((U+du[j])/(omega*r[j]));
	        double alpha = incidence[j]*M_PI/180.0 - phi;
	        double Cl = cl0 + 2*M_PI*alpha;
	        double cla = abs(Cl);
	        double Cd = cd0+cd1*Cl+cd2*Cl*Cl;
	        if (cla > clmax){
	            Cl = clmax*Cl/cla;
	            Cd = cd0 + cla*abs(tan(alpha));
	        }
	        if (Cd > clmax){
	            Cd = clmax;
	        }
	        double gamma = 0.5*V*chord[j]*Cl;
	        gamma = max(gamma,0.0);


	        double Vt = sqrt(pow(U+du[Nstations-1],2) + pow(omega*radius,2));
	        double sinphit = (U+du[Nstations-1])/Vt;
	        double f = Nblades/2.0*(1-r[j]/radius)/sinphit;
	        f = max(f,0.0001);
	        double kappa = 2.0/M_PI*acos(exp(-f));
	        kappa = max(kappa,0.01);
	        kappa = 1; // why?

	        du[j] = sqrt(U*U + Nblades*gamma*omega/(M_PI*kappa))/2.0 - U/2.0;
	        double q = 0.5*V*V;
	        T += (4.0*M_PI*r[j]*du[j]*kappa - Nblades*q*chord[j]*Cd/V)*(U+du[j])*dr;
	        Q += (Nblades*gamma*(U+du[j]) + Nblades*q*chord[j]*Cd*omega*r[j]/V)*r[j]*dr;
	    }

	    double Pout_p = T*U;
	    double Pin_p = omega*Q;
	    eta = Pout_p/Pin_p;
	    CT = T/(pow(RPM/60,2)*pow(2*radius,4));
	    CQ = Q/(pow(RPM/60,2)*pow(2*radius,5));
	    CP = Pin_p/(pow(RPM/60,3)*pow(2*radius,5));
	}
}

double Propeller::getEfficiency(double Jnew){
	if (Jnew != J) setAdvanceRatio(Jnew);
	return eta;
}

double Propeller::getThrustCoefficient(double Jnew){
	if (Jnew != J) setAdvanceRatio(Jnew);
	return CT;
}

double Propeller::getTorqueCoefficient(double Jnew){
	if (Jnew != J) setAdvanceRatio(Jnew);
	return CQ;
}
double Propeller::getPowerCoefficient(double Jnew){
	if (Jnew != J) setAdvanceRatio(Jnew);
	return CP;
}
